#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include "make_ng86.h"

#define NG86_TYPE "char"

int main(int argc, char *argv[])
{
	int i, j, k, l;
	int codon_index;
	struct dsresult value, result[64][64];
	FILE *fh;

	if (argc == 1) {
		fprintf(stderr, "Usage: %s <file>\n", argv[0]);
		return 1;
	}
	for(j = 0; j < 4; j++) {
		for(k = 0; k < 4; k++) {
			for(l = 0; l < 4; l++) {
				codon_index = (j * 4 + k) * 4 + l;
				seq[codon_index][0] = j;
				seq[codon_index][1] = k;
				seq[codon_index][2] = l;
			}
		}
	}

	for(i = 0; i < 64; i++) {
		for(j = i; j < 64; j++) {
			value = calc_ds(3, i, j);
			result[i][j] = value;
			result[j][i] = value;
		}
	}

	fh = fopen(argv[1], "w");

	fprintf(fh, "#ifndef NG86_DS_H_\n");
	fprintf(fh, "#define NG86_DS_H_\n\n");
	fprintf(fh, "/* This file is generated by make_ng86\n");
	fprintf(fh, " *\n");
	fprintf(fh, " * Author: Du Wenfeng <duwf@genomics.org.cn>\n");
	fprintf(fh, " *         modified by Li Heng <lh3@sanger.ac.uk>\n");
	fprintf(fh, " */\n\n");

	fprintf(fh, "const float pot_smut[64] = { /* potential same mut */\n");
	for(i = 0; i < 4; i++) {
		for(j = 0; j < 4; j++) {
			for(k = 0; k < 4; k++) {
				fprintf(fh, "%1.5f", result[0][(i * 4 + j) * 4 + k].s2);
				if (i != 3 || j != 3 || k != 3) fputc(',', fh);
			}
		}
		fprintf(fh, "\n");
	}
	fprintf(fh,"};\n");

	fprintf(fh, "const %s diff_smut[64][64] = { /* same diff */\n", NG86_TYPE);
	for(i = 0; i < 64; i++) {
		fprintf(fh, "{");
		for(j = 0; j < 64; j++) {
			if(result[i][j].ns < 0) result[i][j].ns = -1;
			fprintf(fh, "%d", result[i][j].ns);
			if (j < 63) fputc(',', fh);
		}
		if (i < 63) fprintf(fh,"},\n");
		else fprintf(fh,"}\n");
	}
	fprintf(fh,"};\n");

	fprintf(fh, "const %s diff_nmut[64][64] = { /* not same diff */\n", NG86_TYPE);
	for(i = 0; i < 64; i++) {
		fprintf(fh, "{");
		for(j = 0; j < 64; j++) {
			if(result[i][j].nn < 0) result[i][j].nn = -1;
			fprintf(fh, "%d", result[i][j].nn);
			if (j < 63) fputc(',', fh);
		}
		if (i < 63) fprintf(fh,"},\n");
		else fprintf(fh,"}\n");
	}
	fprintf(fh,"};\n");

	fprintf(fh, "const %s diff_path[64][64] = { /* path */\n", NG86_TYPE);
	for(i = 0; i < 64; i++) {
		fprintf(fh, "{");
		for(j = 0; j < 64; j++) {
			if(result[i][j].np < 0) result[i][j].np = -1;
			fprintf(fh, "%d", result[i][j].np);
			if (j < 63) fputc(',', fh);
		}
		if (i < 63) fprintf(fh,"},\n");
		else fprintf(fh,"}\n");
	}
	fprintf(fh,"};\n");

	fprintf(fh, "#endif\n");

	fclose(fh);
	return 0;
}

int substitution(int codon[], int result[][3][3])
{
	int j, sub_pos, sub_base, na;
	for(sub_pos = 0; sub_pos < 3; sub_pos++) {
		j = 0;
		for(sub_base = 0; sub_base < 4; sub_base++) {
			if(sub_base != codon[sub_pos]) {
				for(na = 0; na < 3; na++) {
					result[sub_pos][j][na] = codon[na];
				}
				result[sub_pos][j][sub_pos] = sub_base;
				j++;
			}
		}
	}
	return 0;
}

struct dsresult calc_ds(int length, int seq_index_0, int seq_index_1)
{
	int s, n, p, i;
	int codon[2][3];

	double s_total_0 = 0, s_total_1 = 0;
	double sd_total, nd_total;

	struct dsresult result;

	sd_total = 0.0;
	nd_total = 0.0;

	i = 0;

	readcodon(codon[0], seq_index_0, i);
	readcodon(codon[1], seq_index_1, i);

	potential_mut(codon[0], &s_total_0);
	potential_mut(codon[1], &s_total_1);

	s = 0;		n = 0;		p = 0;
	mut(codon[0], codon[1], &s, 0, &n, 0, &p);

	result.nn = n;
	result.np = p;
	result.ns = s;
	result.s1 = s_total_0;
	result.s2 = s_total_1;
	return result;
}

int readcodon(int codon[], int seqindex, int pos) 
{
	int k, cvalue;
	for(k = 0; k < 3; k++) {
		codon[k] = seq[seqindex][pos + k];
		if(codon[k] == 4) return 1;
	}
	cvalue = calc_codon(codon);
	return 0;
}

void potential_mut(int codon[], double *s_total)
{
	int k, l, same, nonsense, substcodon[3][3][3];

	substitution(codon, substcodon);
	for(k = 0; k < 3; k++) {
		same = 0;
		nonsense = 0;
		for(l = 0; l < 3; l++) {
			if(aa[calc_codon(substcodon[k][l])] == stopcodon) {
				nonsense++;
			} else {
				if(aa[calc_codon(substcodon[k][l])] == aa[calc_codon(codon)]) same++;
			}
		}
		*s_total += (double)same / (double)(3 - nonsense);
	}
}

int mut(int codon_1[], int codon_2[], int *s_all, int s_path, int *n_all, int n_path, int *p)
{
	int temp[3], codon;
	int i, j;

	if(calc_codon(codon_1) == calc_codon(codon_2)) {
		*s_all += s_path;
		*n_all += n_path;
		(*p)++;
		return 0;
	}
	for(i = 0; i < 3;i++) {
		if(codon_1[i] != codon_2[i]) {
			for(j = 0; j < 3; j++) {
				temp[j] = codon_1[j];
			}
			temp[i] = codon_2[i];
			codon = calc_codon(temp);
			if(aa[codon] != stopcodon) {
				if(aa[calc_codon(codon_1)] == aa[codon]) {
					mut(temp, codon_2, s_all, s_path + 1, n_all, n_path, p);
				} else {
					mut(temp, codon_2, s_all, s_path, n_all, n_path + 1, p);
				}
			}
		}
	}
	return 0;
}

int calc_codon(int codon[])
{
	return ((codon[0] * 4 + codon[1]) * 4 + codon[2]);
}
